Ссылка на практикум
Для выполнения практикума необходимо установить необходимые переменные:
export PATH=${PATH}:/home/preps/golovin/progs/bin
export MOPAC_LICENSE=/home/preps/golovin/progs/bin
- Создадим файл со SMILES-записью порфирина:
%%bash
echo "c1cc2cc3ccc(cc4ccc(cc5ccc(cc1n2)[nH]5)n4)[nH]3 porphyrin" > 1.smi
С помощью программы из openbabel посторим 3D-модель порфирина:
%%bash
obgen 1.smi > 1.mol 2&> /dev/null
%%bash
pymol 1.mol
2 последних водорода (39-ый и 40-ой) оказались лишними.
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys
# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]
### Если вывод в графическое окно тормозит или не нужен, то:
##__main__.pymol_argv = [ 'pymol', '-cp' ]
import pymol
pymol.finish_launching()
from pymol import cmd
### Информацию об ошибках можно смотреть там где запускали ipython notebook
### Будем вставлять файлы изображений
from IPython.display import Image
def Img(name):
cmd.ray()
cmd.png(name)
cmd.delete('all')
cmd.load("1.mol")
cmd.label(selection='id 39-40', expression='" %s" % ID')
cmd.ray()
cmd.png("1.png")
Вот они, торчат из плоскости:
Image(filename='1.png')
Уберём их:
cmd.remove("id 39-40")
cmd.save("1.pdb")
#cmd.delete('all')
#cmd.load("1.pdb")
cmd.ray()
cmd.png('2.png')
Image(filename='2.png')
Готово. Однако он до сих пор не плоский:
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png('2_2.png')
Image(filename='2_2.png')
Теперь переформатируем координаты в mol формате во входной файл для Mopac (c типом параметризации "pm6"):
%%bash
babel -ipdb 1.pdb -omop 1_opt.mop -xk "PM6"
Запустим Mopac:
%%bash
MOPAC2009.exe 1_opt.mop
Переведём результат в pdb:
%%bash
babel -imopout 1_opt.out -opdb 1_opt.pdb
cmd.delete('all')
cmd.load('1_opt.pdb')
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png("1_opt.png")
Он стал плоским!
Image(filename='1_opt.png')
Попробуем провести оптипизацию "AM1":
%%bash
babel -ipdb 1.pdb -omop 2_opt.mop -xk "AM1"
MOPAC2009.exe 2_opt.mop
babel -imopout 2_opt.out -opdb 2_opt.pdb
cmd.delete('all')
cmd.load('2_opt.pdb')
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png("2_opt.png")
Визуально разницы не видно:
Image(filename='2_opt.png')
Наложим изображения друг на друга:
cmd.delete('all')
cmd.load('1_opt.pdb')
cmd.color('tv_red', '1_opt')
cmd.load('2_opt.pdb')
cmd.color('tv_blue', '2_opt')
cmd.ray()
cmd.png("sup_opt.png")
Image(filename='sup_opt.png')
При увеличении видно, что структуры немного расходятся, особенно водороды. Значит в одной из оптимизаций структура получается менее плоской.
- Рассчитаем возбужденные состояния порфирина и спектр поглощения молекулы.
%%bash
cp 1_opt.mop 1_opt_spectr.mop
echo "" >> 1_opt_spectr.mop
echo "cis c.i.=4 meci oldgeo" >> 1_opt_spectr.mop
echo "some description" >> 1_opt_spectr.mop
MOPAC2009.exe 1_opt_spectr.mop
f = open('1_opt_spectr.out', 'r')
lines = f.readlines()[-15:-7]
f.close()
energies = [float(line.split()[1]) for line in lines]
print energies, 'eV'
Получили энергии электронных переходов. Рассчитаем соотвутствующие им длины волн и частоты по формулам: \[\nu (PHz)=\frac{E(eV)}{h(eV \cdot s)}\] \[\lambda(nm)=\frac{h(eV \cdot s)c(m/s)}{E(eV)}\]
wavelengths = [1239.84193/e for e in energies]
frequencies = [e/4.135667516 for e in energies]
print 'frequncy, THz\t wavelength, nm\tenergy, eV'
for i in range(len(energies)):
print '{:^13.0f}\t{:^15.0f}\t{:^10.2f}'.format(frequencies[i]*1000, wavelengths[i], energies[i])
- Определим геометрию молекулы O=C1C=CC(=O)C=C1 (бензохинона):
%%bash
echo 'O=C1C=CC(=O)C=C1 benzoquinone' > quinone.smi
obgen quinone.smi > quinone.mol 2&> /dev/null
cmd.reinitialize()
cmd.load('quinone.mol')
cmd.save("quinone.pdb")
cmd.show('spheres', 'all')
cmd.show('sticks', 'all')
cmd.do("set valence, on")
cmd.do("set_bond stick_radius, 0.14, (all), (all)")
cmd.do("set sphere_scale, 0.25, (all)")
cmd.ray()
cmd.png('quinone_obgen.png')
Image(filename='quinone_obgen.png')
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png('quinone_obgen2.png')
Image(filename='quinone_obgen2.png')
Он плоский. Даже не хочется запускать Mopac..
%%bash
babel -ipdb quinone.pdb -omop quinone_pm6.mop -xk "PM6"
MOPAC2009.exe quinone_pm6.mop
babel -imopout quinone_pm6.out -opdb quinone_pm6.pdb
cmd.delete('all')
cmd.load('quinone_pm6.pdb')
cmd.do("set valence, on")
cmd.ray()
cmd.png('quinone_mopac.png')
Image(filename='quinone_mopac.png')
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png('quinone_mopac2.png')
Image(filename='quinone_mopac2.png')
cmd.reinitialize()
cmd.load('quinone_pm6.pdb')
cmd.color('tv_red', 'quinone_pm6')
cmd.load('quinone.pdb')
cmd.color('tv_blue', 'quinone')
cmd.ray()
cmd.png('quinone_align.png')
Image(filename='quinone_align.png')
cmd.rotate(axis='x', angle=90)
cmd.ray()
cmd.png('quinone_align2.png')
Image(filename='quinone_align2.png')
Почти идентичны. Однако при пристальном изучении можно заметить, что неоптимизированная структура чуть менее плоская.
Повторим всё для дианиона:
f = open('quinone_pm6.mop', 'r')
content = f.read()
f.close()
content = content.replace('PM6', 'PM6 CHARGE=-2')
content = content.replace('O', 'O(-)')
g = open('quinone_pm6_anion.mop', 'w')
g.write(content)
g.close()
%%bash
MOPAC2009.exe quinone_pm6_anion.mop
babel -imopout quinone_pm6_anion.out -opdb quinone_pm6_anion.pdb
cmd.reinitialize()
cmd.load('quinone_pm6_anion.pdb')
cmd.color('tv_green', 'quinone_pm6_anion')
cmd.load('quinone_pm6.pdb')
cmd.color('tv_red', 'quinone_pm6')
cmd.ray()
cmd.png('quinone_align3.png')
Image(filename='quinone_align3.png')
Хмм...
cmd.rotate(axis='x', angle=95)
cmd.ray()
cmd.png('quinone_align4.png')
Image(filename='quinone_align4.png')
Структура плоской быть не перестала. Однако она вытянулась вдоль оси, проходящей по атомам кислорода: сами кислороды разъехались друг от друга (вестимо, из-за взаимодействия зарядов с \(\pi\)-системой), а бензольное кольцо сузилось. Но эти изменения, как мне кажется, не должны сильно превышать изменения от тепловых колебаний.
- Попробуем увидеть переход из тиминового димера в тимины при возбуждении системы:
%%bash
wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/td.pdb" -O td.pdb
cmd.reinitialize()
cmd.load('td.pdb')
cmd.do('set valence, on')
Img('td.png')
Image(filename='td.png')
Оптимизируем структуру с параметрами PM6 при зараде 0:
%%bash
babel -ipdb td.pdb -omop td.mop -xk "PM6"
MOPAC2009.exe td.mop
babel -imopout td.out -opdb td_opt.pdb
cmd.reinitialize()
cmd.load('td_opt.pdb')
cmd.do('set valence, on')
Img('td_opt.png')
Image(filename='td_opt.png')
Оптимизируем результат при заряде +2 (имитируем возбуждённое состояние):
%%bash
babel -ipdb td_opt.pdb -omop td_opt2.mop -xk "PM6"
sed -i 's/PM6/PM6 CHARGE=+2/' td_opt2.mop
MOPAC2009.exe td_opt2.mop
babel -imopout td_opt2.out -opdb td_opt2.pdb
cmd.reinitialize()
cmd.load('td_opt2.pdb')
cmd.do('set valence, on')
Img('td_opt2.png')
Image(filename='td_opt2.png')
%%bash
babel -ipdb td_opt2.pdb -omop td_opt0.mop -xk "PM6"
sed -i 's/PM6 CHARGE=+2/PM6/' td_opt0.mop
MOPAC2009.exe td_opt0.mop
babel -imopout td_opt0.out -opdb td_opt0.pdb
cmd.reinitialize()
cmd.load('td_opt0.pdb')
cmd.do('set valence, on')
cmd.rotate(axis='y', angle='-20')
Img('td_opt0.png')
Image(filename='td_opt0.png')
Как следует из рисунков, при возбуждении системы полного разрушения тиминового димера не произошло. Однако, если я правильно могу интерпретировать результаты, возбуждённое состояние оказалось переходным, из которого система может свалиться как в тиминовые нуклеотиды, так и обратно в димер. А так как система из двух тиминов (-3273.69661 eV) оказалась немного выгодней системы из димера (-3273.58217 eV), то она свалилась именно в тимины.
- Попробуем восстановить ион магния в структуре:
%%bash
wget "http://kodomo.fbb.msu.ru/FBB/year_08/term6/test.pdb" -O atp.pdb
cmd.reinitialize()
cmd.load('atp.pdb')
Img('atp.png')
Image(filename='atp.png')
Добавим к структуре водороды при pH=7:
%%bash
babel -ipdb atp.pdb -omop atp.mop -xk "PM6" -p 7
babel -imop atp.mop -opdb atpH.pdb
cmd.reinitialize()
cmd.load('atpH.pdb')
cmd.set_view ('\
0.850807786, 0.185234621, -0.491741031,\
-0.098345406, 0.975402892, 0.197270557,\
0.516186953, -0.119479574, 0.848098159,\
0.000000000, 0.000000000, -40.559150696,\
-27.376327515, 24.080129623, 6.675564289,\
31.977142334, 49.141159058, -20.000000000' )
Img('atpH.png')
Image(filename='atpH.png')
Добавим к структуре магний:
%%bash
cp atp.mop atpMg.mop
sed -i 's/O -29.24900 1 25.83400 1 5.04800 1/O -29.24900 1 25.83400 1 5.04800 1\nMg -29.15500 1 26.47700 1 5.32550 1/' atpMg.mop
babel -imop atpMg.mop -opdb atpMg.pdb
cmd.reinitialize()
cmd.load('atpMg.pdb')
cmd.label(selection='all', expression='" %s" % ID')
Img('atpMg.png')
Image(filename='atpMg.png')
Запретим движение для всех атомов кроме \(\alpha, \gamma\)-фосфатов, воды и магния:
f = open('atpMg.mop', 'r')
content = f.readlines()
f.close()
s = ''
atomsIDs = [9,41,61,62,12,25,29,33,10,23,27,31,32,14]
for index, line in enumerate(content):
if index<3:
s += line
elif index-2 not in atomsIDs:
temp = line.split()
s += "{:3}{:9} 0 {:8} 0 {:>8} 0\n".format(temp[0], temp[1], temp[3], temp[5])
else:
s += line
g = open('atpMgfreez.mop', 'w')
g.write(s)
g.close()
Сделаем заряд структуры -3:
%%bash
cp atpMgfreez.mop atpMgfreez-without_charge.mop
sed -i 's/PM6/PM6 CHARGE=-3/' atpMgfreez.mop
Запустим оптимизацию:
%%bash
MOPAC2009.exe atpMgfreez.mop
babel -imopout atpMgfreez.out -opdb atpMgfreez.pdb
cmd.reinitialize()
cmd.load('atpMgfreez.pdb')
cmd.distance("id 9","id 41")
cmd.distance("id 9","id 8")
cmd.distance("id 9","id 27")
cmd.distance("id 9","id 29")
cmd.distance("id 9","id 24")
cmd.show('spheres', 'id 9')
cmd.do("set sphere_scale, 0.25, (id 9)")
cmd.set_view ('\
0.850807786, 0.185234621, -0.491741031,\
-0.098345406, 0.975402892, 0.197270557,\
0.516186953, -0.119479574, 0.848098159,\
0.000000000, 0.000000000, -40.559150696,\
-27.376327515, 24.080129623, 6.675564289,\
31.977142334, 49.141159058, -20.000000000' )
Img('atpMgfreez.png')
Image(filename='atpMgfreez.png')
cmd.load('atpMg.pdb')
cmd.color('cyan', 'atpMg')
Img('atpMgalign.png')
Image(filename='atpMgalign.png')
Как видно из последнего рисунка, при оптимизации передвинулись атомы Mg, воды и \(\gamma\)-фосфата. Атомы \(\alpha\)-фосфата остались на месте, не смотря на то, что им было разрешено двигаться.
Ион магния занял положение, при котором он координируется 5-ю связями (с водой, кислородом аспартата и кислородами фосфатов).
А что если Mopac не сообщить информацию о заряде системы:
%%bash
MOPAC2009.exe atpMgfreez-without_charge.mop
babel -imopout atpMgfreez-without_charge.out -opdb atpMgfreez-without_charge.pdb
cmd.reinitialize()
cmd.load('atpMgfreez-without_charge.pdb')
cmd.distance("id 9","id 41")
cmd.distance("id 9","id 8")
cmd.distance("id 9","id 27")
cmd.distance("id 9","id 29")
cmd.distance("id 9","id 24")
cmd.show('spheres', 'id 9')
cmd.do("set sphere_scale, 0.25, (id 9)")
cmd.set_view ('\
0.850807786, 0.185234621, -0.491741031,\
-0.098345406, 0.975402892, 0.197270557,\
0.516186953, -0.119479574, 0.848098159,\
0.000000000, 0.000000000, -40.559150696,\
-27.376327515, 24.080129623, 6.675564289,\
31.977142334, 49.141159058, -20.000000000' )
cmd.load('atpMgfreez.pdb')
cmd.color('cyan', 'atpMgfreez')
Img('atpMgalign2.png')
Image(filename='atpMgalign2.png')
В общем, мягко говоря, результаты изменились. Трифосфатная группа неестественно изогнулась, вода переместилась, кислороды \(\gamma\)-фосфата сблизилсь чересчур близко (т.к. на них нет зарядов и они слабее отталкиваются), однако магний всё равно занял своё законное положение.